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Figure 1: Conjectured phase diagram of QCD as a function of quark chemical potential jJ. and temperature 
T, from Wikipedia. 



1. Introduction 



Just like with water, the form taken by quark matter depends on its temperature and its density, 
or equivalently the chemical potential coupled to the quark number. In fact, one should consider 
distinct, possibly different chemical potentials coupled to the densities of u,d, and s quarks which 
are separately conserved by the strong interactions, giving in total a 4-dimensional parameter space 
with a rich phase diagram. A conjectured two-dimensional (/x, T) section (where the requirements 
of electric neutrality and of beta-equilibrium reduce the three chemical potentials to a single combi- 
nation) is proposed Fig. |l|, taken from a popular reference. The behaviour of QCD in some limiting 
cases (high T or large /x) can be predicted from perturbation theory thanks to asymptotic freedom, 
and the T > 0,/x = properties have been well studied on the lattice. Otherwise, almost all of the 
phase diagram Fig. [T] is based on educated guesses awaiting validation. Putting this phase diagram 
on a firm basis is obviously of fundamental importance. 

The QCD phase diagram follows from the non-perturbative properties of the QCD Lagrangian, 
and can in principle be determined unambiguously by lattice simulations. Unfortunately, standard 
Monte Carlo simulations can only be applied to the /x = vertical axis in Fig. |l]. As is well-known, 
for /I / the simulations are plagued by the "sign problem". The purpose of this review is to 
explain the origin and the nature of the sign problem, and to report on recent progress towards 
circumventing it. I have tried both to start from elementary considerations and to cover some very 
recent, promising developments. This has forced me to skip reviewing some new work presented 
at this Conference, for which I apologize. 
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2. Sign problem 

The sign problem is a necessary evil, unavoidable as soon as one integrates out the fermion 
fields and expresses the partition function in terms of the gauge fields. Analytic integration over 
each fermion species gives a factor det(^ + m + jxyo), where Ip is the massless Dirac operator and 
the last term appears when the chemical potential jJ. is non-zero. Now, Ip satisfies 75-hermiticity: 
75^75 =-^^sothat 

Y5{Ip + m + HYo)Y5=Ip^ + m-HYo = {Ip + m-n*Yoy (2.1) 

Taking the determinant on both sides gives det(^ + m + hyo) = det* (^ + m — /x*7o), which con- 
strains the determinant to be real only if jj. is zero or pure imaginary. In such case, an even number 
of degenerate flavors (same m, same /x) yields a non-negative factor in the integration measure over 
the gauge fields, and standard Monte Carlo techniques apply. The same is true for pairs of flavors 
with opposite, real /x, i.e. isospin chemical potential. 

In the general case, the determinant may be complex, and in fact it must be complex to produce 
the expected physics. This can be seen by considering the free energy of a static color charge 
or anti-charge, respectively related to the expectation value of the Polyakov loop or its adjoint. 
Denoting by dUJ the integration measure which includes the determinant, one sees that 

(Tr Polyakov ) = exp(-}F^) = / Re(Polyakov) x Re(<iGJ)-Im(Polyakov) x Im(<iGJ) (2.2) 
(Tr Polyakov*) = exp(-^F^) = / Re(Polyakov) x Re (<iGT)+Im (Polyakov) x lm{d(n) (2.3) 

Different free energies Fq and Fq, as happens when a chemical potential favors charge over anti- 
charge, can only be obtained if lm{d(n) / 0, i.e. with a complex measure' . 

A corollary of the above statement is that any Monte Carlo ensemble (which is sampled using 
a real non-negative measure) has average baryon number zero (or pure imaginary). So the direct 
sampling of a finite-density ensemble is not possible. Three main approaches have been pursued to- 
wards circumventing this difficulty: reweighting, Taylor expansion and analytic continuation from 
imaginary jj.. I will review them in succession, emphasizing their application to the determination 
of the pseudo-critical temperature Tc{lJ.) and of the QCD critical point. 

3. Reweighting 

3.1 General results 

Let me first illustrate the problem in a toy model. Consider the "partition function" Z(A) = 
J^^ dxexp{—x^ + /Ax). Since Z(A) is real, we can focus on the real part of the integrand, shown 
Fig. ^. While the important values of x are clearly near zero when A = 0, this is no longer true when 
A / 0. Large cancellations take place, and integration far into the tail of the distribution is needed to 
obtain the analytic result Z(A)/Z(0) = exp(— A^/4). The size of the "important" integration region 
is governed by A, not by the width of the A = Gaussian. The situation is similar in QCD, where 



In the SU(2),Nf = 2 case, the square of the determinant remains real positive even when fl ^0. But the baryonic 
chemical potential can be turned into an isospin chemical potential by a redefinition of the quark fields. 
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Figure 2: Toy example of an oscillatory integrand: "every x is important". 



configurations suppressed by i^(exp(— Volume)) must be properly summed over, as stressed most 
forcefully by Splittorff and collaborators [|I]]. Loosely speaking, "every configuration is important", 
and it is not even cleai^ how to sample them. 

In general, one may ask the question: given a real, oscillatory integrand f{x), what is the 
optimal weight g{x) > which should be used for sampling? Uniform sampling is not the answer, 
since it is clearly wasteful to sample regions where |/(x)| is zero or very small. A precise answer 
can be obtained by considering how one forms the expectation value {W)f of a general observable 
W in the desired, target ensemble with partition function Zf = Jdx f{x): 



m . 



{W)f: 



jdxW{x)f{x) _ SdxW{xy-^^g{x) _ {WL 



Jdxf{x) 



Jdx^gix) 



(I). 



(3.1) 



This is the strategy of reweighting: successive measurements of W, obtained from ordinary Monte 
Cai^lo sampling of the auxiliary partition function Zg = Jdx g{x), g{x) > 0, are given a varying, 
oscillatory weight f/g in the ensemble average. The denominator {^)g = Zf/Zg is called the 
"average sign". As we will see shortly, it becomes exponentially small as the volume is increased. 
In addition, the relative eiTor on the average sign propagates to every observable W. Therefore, g{x) 
should be chosen so as to minimize the relative variance of f/g. In the limit where the average 
sign tends to zero, the solution is g{x) = \f{x)\ (up to an ai^bitrary multiplicative constant) |Q]. 
Then, each measurement of the reweighting factor f/g gives ±1, and Zg is often called the "sign 
quenched" ensemble. 

Since the average sign is a ratio of two partition functions, it can be rewritten as 



/ 



-^ = exp 

Zg v^~ ^ 



-^A/(r,A) 



(3.2) 



for a system of volume V at temperature T, where A/ is the free energy density difference be- 
tween the two ensembles, which depends on the temperature and the couplings of the theory (/i 
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Figure 3: Left: Sketch of the QCD pseudo-critical line Tc{pL) (in red), starting from ~ mi^/'i at T = 0, super- 
imposed with the phase transition line (in blue) of the phase-quenched theory (alias isospin chemical poten- 
tial), starting from m-n/l at T = 0. Bose-Einstein condensation of charged pions in the phase-quenched en- 
semble causes a severe sign problem. Right: Comparison of values of the "average phase factor" (exp(2i ) ) , 
measured in lattice simulations and predicted by one-loop chiral perturbation theory [||]. Good agreement 
persists up to T /T^ ^ 0.90. 



for QCD). This makes clear the dependence of the average sign on the volume. To maintain sta 
tistical accuracy, the numl 
exponentially fast with V. 



tistical accuracy, the number of independent measurements must grow as exp(+2^A/(r,A)), i.e.. 



3.2 Reweighting for QCD 

For QCD, the optimal choice of Monte Carlo probability is therefore |Re(det(/x)'^')| ^ |3|]. 
However, this expression cannot be recast, as customary, into a Gaussian integral for further 
stochastic estimation. This causes a considerable ^(V^) computational overhead. A more ap- 
propriate choice is the "phase quenched" ensemble with probability |det(ju)'^/|. Given eq. (^), 
this can be rewritten as det(+/i)^/'^det(— /i)''''/^, corresponding to an isospin chemical potential 
applied to Nf/2 pairs of flavors (/, j). This latter form can be readily recast into a Gaussian integral 
and sampled with the usual Rational Hybrid Monte Carlo. However, the chemical potential is now 
coupled to all qiqj mesons. The lightest among those, the charged pions uy^d (using Nf = 2 nota- 
tion), undergo Bose-Einstein condensation when fi > l~ic{T), with /ic(r = 0) = m-^/l, as illustrated 
Fig. m left. Then, the physics of the auxiliary Monte Carlo ensemble differs qualitatively from that 



of the target baryonic-/^ ensemble, which causes A/ in eq. ( |3.2| ) to become large: the sign problem 
becomes severe. 

Remarkably, the average sign is mostly determined by the physics of pions, for which chiral 
perturbation theory, or even random matrix theory, provides an accurate analytic description pro- 
vided T <^ rriji and fi <C mp/2. A convenient observable to study is the "average phase factor" 

(exp(2/0)) = z(+u'-u) ~ ^ ldetlu!l^ )ll' where Z(+/i,+/i) is the target ensemble with chemical po- 
tential fL, Z{+fL,—fL) is the auxiliary Monte Carlo ensemble with isospin chemical potential, and 
(..)|| is an expectation value with respect to the latter. The observable (exp(2/0)) measures the 
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Figure 4: Left: QCD phase diagram from [0] obtained by combined reweighting in /x and j3 of the /x = 
0,j3 = Pc ensemble (blue dot). Middle: corresponding smallest Lee- Yang zero imaginary part (related to 
the inverse of the specific heat) extrapolated to the thermodynamic limit. Right: full data illustrating the 
insensitivity to ji followed by an abrupt change, courtesy of Z. Fodor. 



effect of changing the chemical potential from +/i to —jJ. for half of the quark flavors, or equiv- 
alently the coiTcsponding fermion boundary conditions in Euclidean time. It is thus ultra-violet 
finite and can be estimated using continuum chiral perturbation theory. Comparison between an- 
alytic and numerical lattice QCD results shows good agreement, even at rather high temperatures 
not far from Tc [^ : see Fig. || right. 

Pion condensation is a consequence of choosing to sample a Monte Carlo ensemble with 
isospin chemical potential. One may wonder if the severe sign problem could not be avoided with 
another choice of Monte Carlo ensemble, e.g., /i = 0. However, it is still because of the phase of 
the reweighting factor that pion condensation for /i > mji/2 does not occur. Additional fluctuations 
in its magnitude only cause the evaluation of the average sign to be even more noisy, and the sign 
problem even more severe. In addition, a pernicious effect occurs: the distribution of the reweight- 
ing factor becomes broader and non-Gaussian, which makes the analysis of the statistical error less 
reliable. Moreover, the necessarily finite Monte Carlo sample may contain zero configurations in 
the region which is most important for the target ensemble. This leads to wrong results, because the 
underestimated statistical error does not reflect the large deviation from the correct answer. Failure 
of reweighting may go undetected, as in the unfortunate "Glasgow method" [^. 

In spite of these difficulties, reweighting has produced a landmark result |^, with a determi- 
nation of the pseudo-critical temperature Tc{lJ.) and of a critical point (see Fig. ^ for QCD with 
physical quark masses, on an A'f = 4 (4 time-slices, a ~ 0.3 fm) lattice. The Monte Carlo ensem- 
ble chosen, (/i = 0,j8 = /3c (/i =0)), was sub-optimal, yet vastly superior to the earlier "Glasgow 
method" which kept j3 fixed [^]. Still, one may question whether the statistical eiTor, which one 
would expect to grow exponentially with /x^, is reliably evaluated. Doubts are fueled by the ob- 
servation of [S] that the critical point is located near the estimated boundary to the pion-condensed 
phase (Fig. ^ left), and that the system appears insensitive to n until very close to the critical point 
(Fig. 0, right). Six years onward, this issue is still not settled, in spite of the authors' own efforts to 
devise a more reliable eiTor estimation [^, and doubts about the determination of the critical point 
Unger on. 

Progress is continuing in the analytic estimation of the average sign. Recent work includes 



the influence of the topological sector at finite /i [12], and the determination of the distribution 
of the phase of the determinant, as well as its correlation with various observables like the baryon 



Finite mu 



Philippe de Forcrand 





Figure 5: Isolines of the average sign in the (lJ.,T) plane for a random matrix model 
sity, j3) plane for Nf = 8 simulations [[1 



and in the (den- 



number []130. Even a random matrix model with a critical point [|10|] gives a description of the 
sign problem roughly consistent with numerical reweighting JTTl]: compare Fig. |5| left and right. 
One should note that pion interactions do not play a major role in the determination of the average 
sign. On the other hand, taking the baryons into account improves the description (and turns out 
to make the sign problem less severe). This can be accomplished by describing the two systems 
(baryonic-/! and isospin-/i) by a non-interacting hadron resonance gas. A fit of lattice simulation 



results to such an ansatz [ 14 1 works well. Interestingly, one can then in turn predict the maximum 
baryon number which can be included in a lattice by reweighting, for an average sign of, say, 0. 1 or 
greater. For practical lattice sizes, this number is ^(10) (and decreases for lighter quarks), which 
is barely sufficient for a statistical treatment. 



4. Taylor expansion 

As we have seen, reweighting is limited to small volumes, and its breakdown is difficult to 
detect. It may be more useful and efficient to try and determine, in the thermodynamic limit, the 
first few Taylor coefficients in the expansion of an arbitrary observable in powers of jx/T about 
/i = 0. In particular, one may consider the pressure P{T,}Jl), since all thermodynamic properties 
can be extracted from its derivatives. Defining AP(r,/x) = P{T,jjl) —P{T,h = 0), one expands 



AP(r,At) 

rp4 



I--(r)(^j 



k=i 



M\ 



2k 



I dii .. 
C2k = (Tr (degree 2k polynomial in I/> , - — ))^=o 



(4.1) 
(4.2) 



where the Taylor coefficients C2k can be expressed as expectation values of traces of matrix poly- 
nomials in the jj. =0 ensemble. The trace of each monomial can then be estimated by the standard 
stochastic averaging over "noise vectors". This strategy looks straightforward, and indeed works 
well at low orders: see Fig. ^ [|l^]. Notice however how the statistical errors grow with the Taylor 
order. Since one must go to higher order as n is increased, to keep the truncation error of the Taylor 
expansion under control, an essential practical question is: how does the work increase with the 
order k ? 

The answer has several parts, and a full complexity analysis has not been carried out yet: 
• The number of terms in the degree 2k polynomial grows approximately as 6^*^ [16]. 



The Taylor coefficient C2k has a finite thermodynamic limit, but the monomials in eq. ( p^ grow 
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Figure 6: First three coefficients in the Taylor expansion of the QCD pressure eq. (4J^) versus T/Tc [ |l5| 



as V^'^, implying large cancellations, which can only be controlled by a large (exponential in k) 
increase with V in the number of noise vectors. This is where the complexity ~ exp(y) of the 
reweighting approach resurfaces. After all, the underlying strategy is an extrapolation in jj. at fixed 
/3, much like the Glasgow method. 

• The distribution of the values whose average yields C2k is less and less Gaussian as k increases. 

• Finite-size effects grow with k, since C2k is analogous to a 2fc-point function. 

These considerations should lead us to expect steady, but slow progress in the expansion to higher 
order. In my opinion, increasing k by one requires increasing computer resources by well over two 
orders of magnitude. The cun^ent state of the art is fcmax = 4 (8^'' order expansion) on an Nt = 6 
lattice 0]. 

In this situation, it is worth exploring alternative methods to obtain the Taylor coefficients. 
These are based on simulations at imaginary chemical potential. 



5. Imaginary jU 



The strategy is simple: perform independent simulations at different values of the imaginary 
chemical potential /i = //i,, fit the results with an ansatz, and analytically continue the ansatz to real 
jU. If the ansatz is polynomial, the fit parameters are the usual Taylor coefficients. Although this 
approach has been used mostly to determine the pseudo-critical temperature Tc{}Jl), it has also been 
applied to the pressure, yielding the same Taylor coefficients C2k as in eq 



^. A recent study [jlj] 

is illustrated Fig. ^. At low temperature, the pressure is best described by a hadron resonance 
gas ansatz. For T > 0.95Tc, this ansatz becomes poor, and a better description is obtained by a 
Taylor expansion, which is sensitive to C(,. Similar observations have been made in Ref. [ |I^ ] on a 



smaller lattice. A technical difference is that Ref. 014[ ] measures only the quark density, i.e. the first 
derivative of the pressure, as a function of imaginary quark and isospin chemical potentials both. 



while Ref. [ ]18[ ] measures all derivatives in /i„ , n^ up to 4th order, but as a function of quark chemical 
potential only. It would make sense to many the two approaches: derivatives up to 4th order are 
easy to compute, and imaginary isospin chemical potential straightforward to implement. Another 
important technical issue should be addressed: how to choose the simulated values of imaginary 
chemical potential and the statistics for each value, so as to maximize the accuracy on a given set of 
Taylor coefficients? Larger values of /I; increase the sensitivity to the desired higher-order terms, 
but also the truncation error in the fitted Taylor polynomial. 
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Figure 7: Imaginary quark density as a function of imaginary quark and isospin chemical potentials, for 
temperatures 0.9Tc (left) and 1.25Tc (right), from [p4[]. The surfaces represent hadron resonance gas and 
polynomial fits, respectively. 



6. Results: Tc{^) 

The pseudo-critical temperature can be determined in a variety of ways. As illustrated Fig. || 
left, all determinations agree for small /i/r. In particular, the curvature ?2 in the Taylor expansion 






'-£'^^(^) 



2k 



(6.1) 



can be determined by the most economical method, which seems to be that of analytic continuation 
from imaginary jj.. While in the past, coarse lattices with Nt =4 and 6 time-slices only have been 
considered, this year has seen remarkable progress, with a first extrapolation to the continuum limit 
for physical quark masses [|^. Using Nt = 4, 6, 8 and 10 lattices and an imaginary-^u method, this 
study confirms earlier indications that the curvature of the pseudo-critical line decreases as a — )• 
iQ] and is small 1 2C ] compared to that of the freeze-out curve, defined as the fireball temperature 
below which inelastic collisions stop taking place and the "chemical", hadronic composition of 
the fireball decay products remains frozen. Since one has a crossover at ^u = 0, Tc depends on the 
observable considered, and so does the curvature. Still, all observables yield a curvature smaller 
than that of the freeze-out curve ?2 ~ 2 [21|. 

This result is of phenomenological importance. As seen Fig. ^ right, a flatter pseudo-critical 
line Tc{lJ.) increases the distance from the putative QCD critical point to the freeze-out curve, 
giving more time for a possible signature of criticality to be washed out as the fireball expands 
before hadronization. Note, however, that the determination of the freeze-out curve is still being 
debated (compare 



11|] with, e.g., [|22|, ^). 
The difficulties of determining subleading coefficients t2k, ^ > 1 in the expansion eq. ( ^ 



has been considered in [g5|], using the imaginary-/^ approach in cases free of sign problem (SU{2) 
and SU{3) with isospin /i) where the analytic continuation to real jj. can be checked against a 
direct determination. For real /i, the pseudo-critical line Tc{lJ.) bends down more and more with 
increasing /x. This indicates that the coefficients t2k are positive. Unfortunately, for imaginary jj. the 
Taylor series is then alternating. Successive contributions largely cancel each other, and the t2k^ 
are poorly determined, as illustrated Fig. |[ Technical issues like the best strategy for choosing 
simulation points and the best choice of fitting ansatz are starting to be explored. 
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Figure 8: (Left) Pseudo-critical temperature determined by various approaches for the same lattice theory 
(4-flavor staggered quarks with mass am ~ 0.05 on an A'r = 4 lattice) [E4|. All approaches agree for IJ./T<1. 
(Right) Effect of a small curvature of the pseudo-critical temperature Tdji): the distance from the putative 
critical point to the freeze-out curve increases, and any signature of the critical point tends to be washed out. 
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Figure 9: Analytic continuation of the pseudo-critical line 71 (/x); for imaginary ji the Taylor series is 



alternating, making the determination of the sub leading Taylor coefficients difficult [ 25 1 . 



7. Results: critical endpoint 



Assume that the phase diagram of QCD features a critical point (jU£, T^), as marked by the star 
in Fig. [l|. This critical point may or may not be chiral. A chiral critical point belongs to the chiral 



critical surface, swept by the /i =0 chiral critical line in the lower left corner of Fig. [10| left as jx is 
turned on. A chiral critical point can be brought to /i=0 by tuning the quark masses, otherwise not. 

A general strategy to locate the QCD critical point, chiral or not, is shown by the first arrow 
Fig. I^ ng/if : one looks for a singularity in /x, keeping the quark masses fixed. All such approaches, 
except for the reweighting study of Fodor and Katz [|^, ^, also keep the temperature fixed. One 
then has to address the delicate question of how to correctly determine the temperature Te- For 



a chiral critical point, an alternative is to stay on the critical surface (second arrow Fig. IC right), 
where the critical temperature is implicitly defined as a function of /i and the quark masses. I 
review these two strategies in succession. 
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Figure 10: {Left) Order of the jj. —0 finite temperature transition as a function of the light and strange quark 
masses. (Right) Two strategies to reach the chiral critical point: (1) at fixed quark masses or (2) along the 
critical surface. 

7.1 Fixed mass, fixed temperature: "effective radius of convergence" 

The determination of the Taylor expansion of the pressure, eq. (p~T|), provides in principle 
a simple way to estimate the location of the QCD critical point: the Taylor expansion will stop 
converging. The pole at {}Xe,Te) in the second derivative d^P/djx^ = Xq of the pressure will 
govern the divergence of the coefficients. In fact, there is another pole at {—}JLe-,Te), so that 
d^P/dn^{H,T = Te) o^ I /{He - M) + 1/(M£ + M) °^ (mI - M^)"'- It follows that 



^ = hm 

Te fn-~ 



C2n{T) 



C2n+2{T) 



at T 



(7.1) 



Thus, some indication about the QCD critical point can perhaps be obtained for free from the first 



few Taylor coefficients. Indeed, this approach has been followed by Karsch and collaborators [ |15| ] 
and by Gavai and Gupta [p^. The latter group has made strong statements, like "We find the radius 
of convergence of the series at various temperatures, and bound the location of the QCD critical 
point to be Te/Tc ~ 0.94 and He/Te < 0.6" (abstract Ref. [17]). Such strong statements force me 
to balance them with strong words of caution. 

The first and obvious point is that eq. (7J_) concerns the « — )• oo limit of the Taylor coefficients. 
From a small number of low-order coefficients, the ratios that one can form are neither a lower 
nor an upper bound of any sort on the radius of convergence. In fact, [|l^] considers the Taylor 
expansion of Xq rather than that of the pressure itself. This leads to an equivalent expression for the 
radius of convergence 



Xq 
f2 



1 d^P °° /U\2«-2 

:£2n(2n-i)c2„(r)(^) 

n=l 



T^dpi"^ 



Te 



lim < 



2n{2n-\)c2n{T) 



{ln + 2){2n+\)c2„+2{T) 



at T 



{12) 
(7.3) 
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But for n = 1, the "effective radius of convergence" differs from that obtained from eq. (7.1) by a 
factor \/6- Thus, the coincidence of estimates from successive small values of n depends in part on 
the choice of observable to expand. Moreover, the QCD critical point is in the universality class of 
the 3d Ising model, with known, non-trivial critical exponents, and the susceptibility Xq does not 
simply diverge like (/i| — /i^)^^ leading to a modification of eq. (L3) for finite n. 



Another difficulty is that the Taylor coefficients vary with the temperature, and so does the ra- 



dius of convergence of the expansion. At high temperature, we know from Roberge and Weiss [ |27[ ] 
that there is a first-order transition line at jx/T = in/3. This singularity at negative jj.^ is consistent 
with the measured Taylor coefficients, which alternate in sign at high temperature starting from C4. 
At low temperature, a first-order transition occurs at real /i, which should cause high-order Taylor 
coefficients to all be positive. This remains true as T increases, until the critical endpoint of the 
first-order line is reached at T = T^. When T rises above Te, the real singularity due to the critical 
point branches into a pair of conjugate poles in the complex /i plane, as shown by Stephanov in 



IgSp. These complex poles move towards the imaginary /i axis as T increases, and approach the 



origin closer than the QCD critical point [|28[]. Therefore, the determination of Te should not be 
based on the minimization of the radius of convergence, but on the sign behaviour of the Taylor 
coefficients. The two groups use different prescriptions. Karsch et al. determine Te as the highest 
temperature at which all measured Taylor coefficients are positive, as appropriate for a singularity 
at real /i. Gavai and Gupta choose for Te the lowest temperature T < T^ for which all measured 



Taylor coefficients are positive (see Fig. 11 of [ 17]). 



Finally, the strongest reason for skepticism, in my opinion, comes from the original reweight- 
ing study of Fodor and Katz |^. The observable they focused on, the imaginary part of the Lee- 
Yang zero closest to the real axis, extrapolated to the thermodynamic limit, is closely related to the 
inverse of the maximum value of the specific heat. In the small-/x region where reweighting can 
be trusted, the specific heat shows complete insensitivity to jU. It would take a high-order Taylor 
expansion about /i = of the curve Fig. ^ right to capture the zero signaling the QCD critical 
point. Note, as further confirmation of this lack of sensitivity to a putative critical point, that all 
susceptibilities measured in the determination of the pseudo-critical line Tc{lJ.) [|^] reported Sec. ^ 
appear to slightly decrease as /i is turned on. 

7.2 The curvature of the critical surface 

Considering the difficulties and ambiguities of determining the convergence radius of the Tay- 
lor expansion, it is worth pursuing a complementary, broader strategy: the determination of the 
chiral critical surface Fig. [T^ right, following the second arrow. I have been pursuing this approach 
for several years with Owe Philipsen [ ^ , ^ , 31, 32, 33] . We have now reached conclusive results 



on coarse Nt = 4 lattices, using standard staggered fermions and setting aside possible issues with 
taking fractional powers of the fermion determinant. 

The first step was the determination of the /i = critical line in the (m„rf,mj) quark mass 



plane, shown Fig. 1 1 left. No surprises were encountered there. We confirmed that the physical 
point lies in the region where a finite-temperature crossover takes place, not a phase transition. 
And our results were consistent with the tricritical scaling implied by a tricritical point at m„ ^ = 0, 



nis ~ 500 MeV (solid blue curve in Fig. 1 1 left). 
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Figure 11: Left: Chiral critical line at /i = in the (»i„^/,mj) quark mass plane [31 1. Compare with 
expectations Fig. n^/e/f. The two red arrows mark the points at which the curvature of the critical surface 
was measured. Center. Leading effect of a negative curvature. Note the curvature of the deconfinement 
critical surface for heavy quarks OM]. Right: Possible back-bending due to higher-order terms [p5|]. 



The next step has been to measure the variation of the critical quark masses with chemical 
potential, expressed in Taylor series form as 



m,{0) 



f 1^ \ 



i + E'^i^ 



2k 



k=l 



\71t) 



(7.4) 



This endeavour has been pursued at the two mass points marked by red arrows Fig. |TT] left, cor- 
responding to three degenerate flavors and to a strange quark with physical mass, respectively. In 
the Nf = 3 case, we have compared the effect of an infinitesimal imaginary jj., which yields the 
Taylor coefficients ct directly, and that of several finite imaginary /I's with results fitted by a poly- 
nomial, finding good consistency. We have also compared two spatial volumes, 8^ and 12^, finding 



excellent agreement. Our final result [32] is 



-M., .3.3,3) (ii)^-47(20)(iL^ 



m,(0) 



kTJ 



(7.5) 



Since we identify the critical mass nic as that which gives the Binder cumulant B4 = }L ^Vm of 
the quark condensate yy its critical Ising value 1.604.., our observable is built from 4''' derivatives 
of the pressure. Extracting C2 as above then requires the measurement of 8"* derivatives of the pres- 
sure. As indicated Sec. 0, this is the current state of the art. To obtain such results, we accumulated 
about 25 million RHMC trajectories, and started to use the EGEE computing Grid [p6|]. 

In the Nf = 2 + 1 non-degenerate case with physical strange quark mass, we had to tune mu,d 
to values smaller than in nature in order to turn the finite-temperature crossover into a second- 
order phase transition. This forced us to use a larger, 16^, volume and use the computing Grid 



extensively. Our final result [33], based on about 1.5 million trajectories, is 



mg-"(M) 
m^''"(0) 



1-39(8) 



nTJ 



(7.6) 



Thus, in both cases we find that the region of quark masses where a first-order transition takes 
place shrinks as n is turned on, just like in the case of heavy quarks where the sign problem is mild 
and the critical surface can be determined by brute force reweighting [|^, 37]. The chiral critical 
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Figure 12: Left: estimate of the region of first-order transition in the quark mass plane, in physical units, 
from improved (stout-smeared) actions on A', =4 and 6 lattices [Ea]. Right: current status of the measure- 
ment of the curvature of the critical surface on A', =6 lattices. The (^{{p./T)'^) Taylor coefficient is obtained 
from the intercept, and the ^((/x/T)^) coefficient from the slope of the fit. A negative curvature is favored 
as for Nt = 4, but the results are not conclusive yet. 



surface then bends as shown Fig. 1 1 center. As indicated in the figure, this rules out, on A'f = 4 
lattices, the presence of a chiral critical point, at least for ^/T<&{\) where one would expect the 
truncation error of the Taylor expansion to be small. 

This conclusion must be accompanied by important cautionary remarks. First, one may be 
worried about the convergence of the Taylor expansion. Already fox jx/T = 1, the last terms of 
eqs. ( [T^ ) and ( f/^ are dominant. Even though the sign of the next-order contribution has been 
estimated and reinforces the shrinking of the first-order region, it is clear that higher-order terms 
may quickly produce a "back-bending" of the critical surface as in Fig. [IT] right. Indeed, this must 
happen if one wishes to reconcile our result with that of Fodor and Katz [^], whose critical point 
Fig. ^ lies at jU/r ~ 0.7 for mass parameters similar to eq. ([7^). Similarly, there is no inconsistency 
between our conclusion and that of Ejiri [38], who finds a critical point at pi/T ~ 2.4. Such a 
back-bending can also be explained by model calculations. In the Nf = 2 case, it is known that a 
restoration of the U{1)a symmetry favors a first-order transition [ p9| , pO| ] . Similarly for A^y = 2 + 1 , a 



back-bending surface can be produced in an NJL-type model, by adding a 't Hooft ?7(1)a -breaking 
term (det^,(l — /s)^^ +h.c.), with a /i -dependent strength [Q ^. A similar pattern can even be 



obtained in a linear sigma model including thermal fluctuations []41[]. 

The second issue is the systematic error caused by the rather coarse lattice spacing: Nt = 4 
implies a ~ 0.3 fm. Even at /i = and on Nt = 8 lattices, discretization errors are presumably 
the cause of the (^(15%) discrepancy between estimates of Tc by Karsch et al. |43] and by Fodor 



et al. []44|]. Here, the critical surface which we study is sensitive to deviations from the Stefan- 
Boltzmann law at high temperature, and even more so to the violation of taste symmetry which 
strongly affects the thermodynamics of the 16 "pions". It should come as no surprise - a posteriori 
- that the critical line in the quark mass plane, using physical coordinates like niji/Tc, seems to move 



Interestingly, back-bending does not necessarily produce a critical point. The critical temperature is determined 
implicitly as one moves on the critical surface, and may decrease to zero, thus terminating the critical surface, before the 
quark masses have reached their physical values |p3]. 
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by ^(100%) as Nt increases from 4 to oo. Early indications came from improving the discretization 
of the Dirac operator at fixed Nt =4: the critical pion mass seemed to decrease by a factor r^ 4 [^] . 
More recently, an even more dramatic shrinking of the first-order region has been reported in [p^], 
estimating from stout-smeared Nt = 4 and 6 simulations that the first-order region is limited to 



quark masses about 10 times smaller than physical (see Fig. [12| left). A less dramatic, but similar 
effect has been seen in [Q] when comparing niji/Tc for the 3-flavor theory on lattices with Nt =4 
{niji/Tc = 1.680(4)) and 6 {mji/Tc = 0.954(12)) time-slices. Clearly, finer lattices tend to make 
the finite-temperature transition smoother. This can in fact be seen in an NJL model by truncating 
the sum over Matsubara frequencies to the first Nt values p6[]. This effect by itself makes it more 
difficult to imagine a critical point at small /i : the physical point lies further from the critical line, 
and the critical surface would have to bend very strongly towards larger quark masses. 

Moreover, one may hope, because the chemical potential does not affect the UV physics, that 
while the critical surface will move significantly towards the origin m„ ,;/ = m^ = as a — )• 0, its 
curvature will vary less. This pious wish needs to be confirmed by numerical simulations, which 
have proved to be challenging. Using the same method as for Nf = 4, we have been simulating an 
18^ X 6 lattice, with 3 degenerate quark flavors at the /i = critical mass. After over two years of 
simulation and a half-million units of Molecular Dynamics time, the current results Fig. |^ right 
are still too noisy to conclude. The curvature of the critical surface is obtained from the /x = 
intercept of the data: a negative value is favored, as for Nt = 4, but the errors are large and the 
(correlated) data do not exclude a positive value. Moreover, a linear fit, including a (/i/r)^ term, 
is preferred over a constant fit. The sign of the {n/T)^ coefficient then reinforces the shrinking of 



the first-order region as in eq. (7.5). However, its magnitude (which is very poorly determined at 
the moment) could make this term dominant over the leading {jjl/tY term as soon as }Jl/T>0A5. 
If this turns out to be the case, then our truncated Taylor expansion can only be trusted up to such 
values. Going to larger pi would require determining higher-order Taylor coefficients, a daunting 
task. 

8. Left out 

I have not covered several recent developments in the numerical study of finite-density QCD. 
• The numerical study of the canonical ensemble, with a fixed number of baryons, has been ex- 



tended from staggered fermions p7| , ^ ^] to Wilson fermions, with several technical refine- 
ments [^, 50, 51]. This is computationally challenging, since the fermion determinant must be 
computed exactly, at a cost proportional to the cube of the matrix size, so that the work is increased 
64-fold on a given lattice size. On a 6^ x 4 lattice, the phase diagram appears to be qualitatively 
different for Nf = 2,3 and 4 flavors, suggesting a possible critical point for Nf = 3. These results 



are presented by Anyi Li in a plenary talk [|52|]. 

• Instead of aiming at a finite density of baryons, one may study T = few-body physics: by 
measuring the interactions between two or three baryons, one can constrain the couplings of an 
effective theory describing nuclear matter. This is a very important and active research direction. 



reviewed last year [53], with two large-scale efforts underway [54, 55]. The sign problem appears 



in the form of a signal-to-noise ratio for baryon correlators degrading exponentially fast in Eu- 



clidean time [56]. But since one works at T = 0, a toolbox of variational trial states can be used 
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to try and isolate the groundstate before the signal dies away [p7|]. The state of the art is reviewed 



by Will Detmold in a plenary talk []58|]. Spectacular results have already been obtained for multi- 



meson states, corresponding to an isospin chemical potential []59|]. 

• Other interesting developments concern non-QCD theories, e.g. 2-color "QCD" [ |60| ] where a va- 
riety of regimes may appear at low temperature, as the chemical potential is increased, or (1 + l)(3f 
Gross-Neveu and NJL models whose phase diagram can be determined analytically [ ]6li |62| ], in- 
cluding inhomogeneous, crystalline phases. These features may well be present in QCD also. 

• An important approach to tackling the sign problem is the "density of states" approach, where the 
partition function is expressed as a one-dimensional integral Z = Jdx p{x), and the computer effort 
can be concentrated on values of ;c where p(x) is more noisy [|63|]. A large-scale effort, taking for ;c 
the gluon action, gave hints of a triple point in the (/x, T) plane [^]. More recently, this approach 
has been espoused, in part, by Ejiri in [^]. 

• Several contributions to this Conference did not fit in the categories above, and I omitted them 
from this subjective review, with apologies. 

9. Prospects 

I have painted a rather bleak landscape, where systematic errors are very large and one should 
only expect slow progress, helped with massive amounts of computer time. Is there nothing more 
exciting in sight? Yes, definitely. Let me emphasize two directions where I consider enthusiasm to 
be justified. Amusingly, they both represent a revival of topics which were "hot" twenty years ago. 

9.1 Worldline formalism and strong coupling limit 

As we have seen Sec. ^, as soon as one integrates out the fermions, one must encounter a 
sign problem with the resulting determinant at finite density. This suggests changing the order of 
integration, and performing at least a partial integration over the gauge fields first [|65|]. Integrating 
out the gauge links seems hopeless, since the plaquette term in the action introduces a complicated 
4-link interaction. Still, this may be simpler than dealing with the sign problem. To start with, 
one can consider the strong-coupling limit jSgauge = 0, where the plaquette term in the action drops 
out. Then, the link integration factorizes into a product of 1-link integrals, which can be performed 
analytically. Only color singlets survive, made of quark fields. The Grassmann integration pro- 
duces only a handful of terms, which represent the hopping of color singlets from site to site. In 
other words, the partition function has been reexpressed as a sum over configurations of loops, 
representing the worldlines of hadronic color singlets. 

At this stage, this reformulation is relatively simple, since it only involves discrete variables, 
and physically appealing. But it is not clear that simulations will be any easier: The fixed number 
of underlying Grassmann variables at each site generates a constraint on the loop configuration 
(the number of qq mesons connected to each site is fixed), which makes local Monte Carlo up- 
dates hopelessly inefficient. And a severe sign problem follows from the fermionic nature of the 
baryons (for an odd number of colors): baryon loops are oriented, and their weight flips sign with 
their orientation. Fortunately, this sign problem can be solved by partial resummation at /i = 0, 
and then remains very mild at ^Li / [|6^]. And a recent Monte Carlo algorithm, the "worm" algo- 
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rithm []67[], is particularly efficient at updating globally such discrete systems with constraints [pq]. 
This efficiency does not even degrade in the chiral limit. 

All this was understood and tried around 1990 [^, ^, and abandoned before the discov- 
ery of the worm algorithm. This crucial algorithmic progress has enabled a complete numerical 
determination of the (/x,r) phase diagram for the strong coupling limit of QCD with staggered 
fermions ||70|], superseding ancient ^n\\ and recent [ [72| ] analytical approximations, and even a pre- 
cise study of nuclear interactions and of "nuclear matter" in this lattice model, all with "tabletop" 
computer resources. 

I am insisting on one particular project because of my personal involvement. But the same 
worldline approach can be applied to other lattice fermion models [^, 74], and to bosonic models 



normally afflicted by a sign problem, after transforming to dual variables [[7^]. It may even be useful 
to "fermionize" a bosonic system, treating the bosons as fermion composites. Recall that this yields 
a simple derivation of Onsager's solution to the 2d Ising model ^Wj. This "new computational 
approach" was reviewed by its main proponent at last year's Conference J/T]]. 

Can all lattice models, in particular QCD at weak coupling, be formulated and efficiently 
simulated as a gas of loops? One technical ingredient for success is the worm algorithm. Its limits 
of applicability are not clear yet. It appears to work well for at least one model, the 2d CP^^ ^ spin 
model [J78|], where cluster algorithms are known to fail. If the random worms can be generalized to 
random surfaces, then one could apply the same treatment to the Yang-Mills part of the action and 
simulate QCD at weak coupling, as a gas of quark loops forming the boundary of gauge surfaces. 
Unfortunately, this step requires a duality transformation, which for a non-Abelian theory gives 



negative weights and/or non-local couplings [79, 



A somewhat less ambitious strategy consists of designing fermion-based, sign-problem free 
actions, with symmetries which ensure a desired effective low energy limit. Such actions can be 
efficiently simulated, even in the massless limit, to address precise questions about the effective 



low energy theory. A variety of models with chiral symmetry can be realized [|8ip. What is still 
missing is a non-Abelian gauge symmetry. 

9.2 Complex Langevin 

The idea of stochastic quantization is to introduce a Langevin evolution for a field in a 
fictitious time T, obeying 

where Tj is a Brownian noise. When the action S[<p] is real, one can prove the existence of an as- 
sociated Fokker-Planck equation and its convergence to the fixed-point distribution oc exp(— 5'[0]), 
so that all observables satisfy 

{Wmr = ^J&^cxp{-SmW[^] (9.2) 

where (..)t is an average over the fictitious Langevin time T. 

What happens when S[^] is complex? The drift force — ^ becomes complex, so that each 
component of the field <p will become complex under the evolution eq. ( |9. 1[ ). One can then com- 
plexify <p into (0^ + i<p'), evolve with the complexified version of eq. (pH]), and monitor the time- 
average {W [0^ + /0^])t. There still is an associated Fokker-Planck equation, but almost nothing 
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is known about stationary solutions. Nevertheless, it turns out that {W [(j)'^ + i0^] )t has the desired 
value 2 1 ^^ exp(— 5 [(/>])VK [^] - sometimes. Some other times, the {^^ + /(/>') distribution in the 
complex plane does not converge to a T-stationary distribution, because (^'^ + i^^) runs to infinity. 
And some other times, convergence is achieved to a distribution giving wrong expectation values 
(W [0^ + /0^] )x- An understanding of sufficient conditions (besides taking a real action) to avoid 



runaways and convergence to the wrong answer was not reached in the 1980s [ ]82| , |83| , |84| , |85| ]. 
Activity on this topic died away. 

It was noticed only recently that a drastic reduction in the discrete stepsize of the complex 
Langevin evolution sufficed to almost eliminate runaways. With an adaptive stepsize [|^] the 
runaways disappeared completely. This opened a new playground where to study convergence. 
Surprisingly, correct answers have been obtained in many cases: toy models with one gauge link 
matrix; Ad complex ^^ theory with a chemical potential, where T w Bose-Einstein condensation 
is reproduced for \i > }Xc [p7|l, as well as the lack of /x -induced effect for fX < fXc (known as the 



Silver Blaze problem" []88[]); and even QCD with chemical potential, in the heavy-dense limit. 



where results are consistent with those of reweighting 089[]. These extraordinary successes should 
be balanced against a short list of failures: the noise 77 in eq. (^?T|), which a priori only needs to 
satisfy (t7(t)|77(t')) = 25(t — t') and can be complex, must in fact be kept real for convergence; 
and real-time quantum evolution still seems to be out of reach ||90|]. Clearly, the approach looks 
promising and its limits must be further explored and understood. 



To get the flavor of the magic at work here, let us consider the simple example of Sec. ^ 
Z(A) = J^2 dxex^{—x^ + iXx). When A = 0, the real Langevin evolution is dx/dx = —Ix + rj. 
When A / 0, the drift force becomes complex and one needs to complexify x into [x + iy). The 
corresponding Langevin evolution becomes 

— {x + iy) = -2{x + iy) + /A + T] (9.3) 

dT 

In this simple case, this equation can be solved analytically. Since Z(A) is a Gaussian integral, the 
stationary distribution of {x,y) is a nice Gaussian shown Fig. |T3[ centered at the complex saddle 
point {x = 0,y = /A/2). It is perhaps not intuitive how expectation values {W{x)) are recovered: 
one must analytically continue W{x) to W{x + iy) to obtain {W{x + iy))T: = {W{x))z- The original, 
negative contributions of some values of x are then reconstructed when one integrates VK (x + iy) 
over y. Note that the j-width of the Gaussian depends on the variance of the imaginary part of 
the Langevin noise: for a real noise, the 3'-width shrinks to zero. But, in this case at least, correct 
answers are obtained for any complex Langevin noise satisfying (t](t)|t](t')) = 25(t — t'). 

This toy example suggests that an analysis of the saddle points of the classical action can be 



fruitful. This is the starting point of a loop-like expansion considered in ||91|], which may shed 
light on the convergence properties of complex Langevin: The Langevin noise distorts the classical 
Gaussian distribution around the saddle points, hopefully in a weak manner. Indeed, one may guess 
that systems with one complex saddle point can perhaps be safely studied by complex Langevin, 
while competition between saddle points (as at a first-order transition) may present a challenge. 
The "safe" category might include the effect of a vacuum angle, as studied in a saddle-point 



formulation in [92]. Gauge theories are more likely to be in the "unsafe" category, due to flat 



directions in the action, which correspond to gauge transformations and extend to complex infinity 
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Oscillatory weight(x) 
Positive weight(x,y) 




Figure 13: Complex Langevin for 1 degree of freedom x — > (x + iy): the oscillatory x-distribution of Fig. GJ 
shown in red and green, becomes a smooth positive Gaussian in the {x,y) plane, centered at the complex 
saddle point /A/2. All moments of the original, oscillatory x-distribution are equal to moments of (x + iy) 
with respect to the Gaussian, positive {x,y) distribution. 



because gauge links are complexified from SU{N) to SL{N,C). Progress is being made towai^ds 
understanding necessary conditions for convergence [p3|]. Finally, I note the intriguing, perhaps 



fruitful analogy between complex Langevin and /T-symmetric quantum mechanics [94] and its 



cousin, complexified classical mechanics []95|]. 



10. Conclusions 

Determining the phase diagram of QCD as a function of temperature and chemical potential is 
an important fundamental goal. Even if clear answers are not available yet, and if the progress in 
our knowledge is likely to be slow, it is definitely worth pursuing the present efforts. Finite-density 
lattice QCD is not just a temporarily fashionable topic: it justifies a sustained research program. 

Three numerical approaches give reliable, consistent results provided the chemical potential 
is small enough: reweighting, Taylor expansion and analytic continuation from imaginary /i. This 
allows crosschecks in the region where the reliability of these methods becomes doubtful. Con- 
fidence in finite-density simulations can be established, so that comparison between QCD and 
models can become reliable and fruitful. For instance, the well-established phase diagram of QCD 
at imaginary jj. is already a useful testing ground for effective models. 

Similarly, the curvature of the pseudo-critical temperature TdlJ.) at /i = is almost under nu- 
merical control, and provides useful phenomenological information. The situation is more delicate 
regai^ding the existence and location of a QCD critical point, which requires venturing to non-zero 
/I. Nevertheless, analytic knowledge about the severity of the sign problem can give us, before 
we start the computation, reliable information about the {pi,T) region which can be explored and 
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the baryon density which can be reached. And it now seems clear that the coarseness of the usual 
A^, = 4 or 6 lattices is the major source of systematic error. Thus, with necessary massive increases 
in computer resources, we can expect slow but steady progress towards the final, continuum limit 
answer. 

This review suggests two possibilities for breakthroughs, or at least for rapid development: 
reversing the order of integration and integrating over the gauge links first; and dealing with the 
sign problem using complex Langevin. Whatever the final scope of these two strategies, one can 
already predict that they will contribute to increasing our knowledge of finite-density QCD at least 
in some limiting regimes of parameters. 

I would like to end by citing Confucius, who knew all about the importance of scholarly 
research: "Real knowledge is to know the extent of one's ignorance". 
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